home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / Regression.rexx < prev    next >
Encoding:
OS/2 REXX Batch file  |  1998-10-07  |  14.6 KB  |  791 lines
  1. /* Regression Statistics */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Math Library needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30)
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19.   /* add to library list */
  20. signal off syntax
  21.  
  22. /* Start Main Routine */
  23. if loadflag=1 then 'Load()'
  24. 'ActivateWindow()'
  25. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  26. colon=pos(":",range)
  27. if colon=0 then do
  28.     'Message "Please select a range before executing this script"'
  29.     'DEFPUBSCREEN("Workbench")'
  30.     exit
  31. end
  32.  
  33. /* Find cell references and cell, column numbers */
  34. start_cell=substr(range,1,colon-1)
  35. end_cell=substr(range,colon+1)
  36. start_row=cellrow(start_cell)
  37. end_row=cellrow(end_cell)
  38. start_col=cellcol(start_cell)
  39. end_col=cellcol(end_cell)
  40. NRows=end_row-start_row+1
  41. NCols=end_col-start_col+1
  42. if NCols>2 then do
  43.     'Message "Only Two Columns allowed!"'
  44.     'DEFPUBSCREEN("Workbench")'
  45.     exit
  46. end
  47.  
  48. /* Ensure we have X and Y in the right columns */
  49. in=2 /* Independent variable x in second column */
  50. de=1 /* Dependent variable y in first column */
  51. an=rtgetstring(,"Assume Dependent Variable is First Column?: Y or N","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  52. if an="" | rtresult=0 then do
  53.     'Message "Aborting!"'
  54.     'DEFPUBSCREEN("Workbench")'
  55.     exit
  56. end
  57. an=left(an,1)
  58. an=Upper(an)
  59. if an="N" then do 
  60.     in=1
  61.     de=2
  62.     end
  63.  
  64. /* Get cell reference for output range */
  65. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCalc"')*/
  66. if out_cell="" then do
  67.     'DEFPUBSCREEN("Workbench")'
  68.     exit
  69. end
  70. if length(out_cell)<2 | datatype(left(out_cell,1),'n') then do
  71.     'Message "Invalid cell reference"'
  72.     'DEFPUBSCREEN("Workbench")'
  73.     exit
  74. end
  75. /* Suppress Screen Redraw to Speed Things Up */
  76. 'Refresh 0'
  77.  
  78. /* Open a small output window on tcalc screen*/
  79. fo=0
  80. CR='0a'x
  81. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  82. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  83.      call writeln(6Info, DisplayMsg)
  84.     fo=1
  85. end
  86. else do
  87.     'Message "TCALC Screen not available for Progress messages"'
  88. end
  89. CALL DELAY(150)
  90.  
  91. /* Get cell references for top cell in each column */
  92. 'SelectCell' start_cell
  93. do col=start_col to end_col
  94.     'GetCursorPos'
  95.     top_cell.col=result
  96.     'Column 1'
  97. end
  98. if fo then call writech(6Info,"Progress...10 ")
  99.  
  100. /* Get labels for later use on output */
  101. 'SelectCell' start_cell
  102. 'GetValue'
  103. testlabel=result
  104. testlabel=strip(testlabel)
  105. if datatype(testlabel,'n')=1 then do
  106.     labelflag=0
  107.     do x=1 to NCols
  108.     title.x="Column "||x
  109.     end
  110. end
  111. else do
  112.     labelflag=1
  113.     NRows=NRows-1
  114.     do x=1 to NCols
  115.     'GetValue'
  116.     str=result
  117.     title.x=translate(strip(str),"_"," ")
  118.     'Column 1'
  119.     end
  120. end
  121. if fo then call writech(6Info,"20 ")
  122.  
  123. /* Get data from cell range */
  124. col=start_col
  125. lav=0
  126. tot=0
  127. count.=0
  128. total.=0
  129. do x=1 to NCols
  130.     'SelectCell' top_cell.col
  131.     if labelflag=1 then 'CursorDown 1'
  132.     do y=1 to NRows
  133.         'GetValue'
  134.         valtest=result
  135.         if datatype(valtest)='NUM' then do
  136.             'GetValue'
  137.             val=result
  138.             data.x.y=val
  139.             tot=tot+val
  140.             total.x=tot
  141.             count.x=1+count.x
  142.         end
  143.         'CursorDown 1'
  144.     end
  145.     col=col+1
  146.     tot=0
  147.     lav=0
  148.     val=0
  149. end
  150. if count.1 ~= count.2 then do
  151.     'Message "There must be equal numbers of values in both columns"'
  152.     'DEFPUBSCREEN("Workbench")'
  153.     exit
  154. end
  155. if fo then call writech(6Info,"30 ")
  156.  
  157. /* Calculate Means */
  158. mean.=0
  159. do x=1 to NCols
  160.     mean.x=total.x/count.x
  161. end
  162.  
  163. /* Calculate Standard deviation and Variance */
  164. dat=0
  165. meenx=0
  166. tot.=0 /* Array holding sums of columns */
  167. totsq.=0 /* Array holding square of column sums */
  168. XY=0 /* Var holding X times Y */
  169. diff.=0 /* Array holding difference between value and column mean */
  170. sum.=0 /* Array holding sum of value minus mean of value squared */
  171. sd.=0 /* Standard deviation array */
  172. var.=0 /* Variance array */
  173. do x=1 to NCols
  174.     sum.x=0
  175.     meenx=mean.x
  176.     do y =1 to count.x
  177.     dat=data.x.y
  178.     tot.x=dat+tot.x
  179.     totsq.x=(totsq.x)+dat**2
  180.     diff.x.y=dat-meenx
  181.     sum.x=(dat-meenx)**2+(sum.x)
  182.     end
  183.     N=(count.x)-1
  184.     var.x=(sum.x)/N
  185.     sd.x=sqrt(var.x)
  186. end
  187. Do y=1 to NRows
  188.     XY=XY+((data.1.y)*(data.2.y))
  189. end
  190. say "XY="XY
  191. if fo then call writech(6Info,"40 ")
  192. N=count.1
  193. summulti=0 /* Sum of product of diffs */
  194. diffMult.=0 /* Array holding values of diffx times diffy */
  195. do y=1 to N
  196.     diffMult.y=(diff.de.y)*(diff.in.y)
  197.     summulti=summulti+(diffMult.y)
  198. end
  199.  
  200. /* Calculate a & b and Predicted Y */
  201. b=0
  202. a=0
  203. sumdifpa=0
  204. Pr.=0 /* Array holding values of predicted values */
  205. difpa.=0 /* Array holding differences between observed and predicted values */
  206. b=summulti/(sum.in)
  207. a=(mean.de)-b*(mean.in)
  208. do y=1 to N
  209.     Pr.y=a+b*(data.in.y)
  210.     difpa.y=(data.de.y)-(Pr.y)
  211.     sumdifpa=sumdifpa+(difpa.y)**2
  212. end
  213. if fo then call writech(6Info,"50 ")
  214.  
  215. /* Calculate standard Deviation */
  216. stdev=0
  217. t=0
  218. stdev=sqrt((sumdifpa/((N-2)*(sum.in))))
  219. t=b/stdev
  220. df=N-2 /* Degrees of Freedom */
  221.  
  222. /* Calculate Pearson r */
  223. rs=0
  224. pr=0
  225. rs=((summulti)**2)/((sum.in)*(sum.de))
  226. pr=sqrt(rs)
  227.  
  228. /* Calculate ANOVA */
  229. SSRegTop=(XY-(((tot.1)*(tot.2))/NRows))**2
  230. SSRegBot=(totsq.1)-(((tot.1)**2)/NRows)
  231. SSReg=SSRegTop/SSRegBot
  232. SSTot=(totsq.2)-((tot.2)**2/NRows)
  233. SSRes=SSTot-SSReg
  234. MeanSqReg=SSReg
  235. MeanSqRes=SSRes/df
  236. FR=MeanSqReg/MeanSqRes
  237. DFB=NCols-1
  238. DFW=NRows-NCols
  239. if fo then call writech(6Info,"60 ")
  240.  
  241. /* Calculate Probability */
  242.  
  243. AC=.0001
  244. EC=.005
  245. EC2=EC+EC
  246. PF=PROBABILITY(FR,DFB,DFW)
  247. PF=1-PF
  248.  
  249. if fo then call writech(6Info,"70 ")
  250. /* P=.95 OR .99 */
  251. FCRIT1=FCRITICAL(.95,DFB,DFW)
  252. FCRIT2=FCRITICAL(.99,DFB,DFW)
  253.  
  254. /* Calculate T distribution function */
  255.  
  256. P=PROBT(t,df)
  257. if t>=0 then P=1-P
  258. PT=P*2
  259. /* Calculate T Critical */
  260. if fo then call writech(6Info,"80 ")
  261. TCRIT1=TCRITICAL(.95,df)
  262. TCRIT2=TCRITICAL(.99,df)
  263. TCRIT3=TCRITICAL(.975,df)
  264. if fo then call writech(6Info,"90 ")
  265. TCRIT4=TCRITICAL(.995,df)
  266. if fo then do
  267.     call writeln(6Info,"100")
  268.     call writeln(6Info,"Writing output to window...")
  269. end
  270. /* Output */
  271. PV=0
  272. 'SelectCell' out_cell
  273. 'ColumnWidth 10'
  274. 'Put "Least Squares Regression"'
  275. 'CursorDown 1'
  276. 'Put "Predicted Values"'
  277. 'CursorDown 1'
  278. do y= 1 to N
  279.     'Put' format(Pr.y,,4)
  280.     'CursorDown 1'
  281. end
  282. 'SelectCell' out_cell
  283. 'CursorDown 1'
  284. 'Column 2'
  285. 'ColumnWidth 25'
  286. 'Put "Regression Statistics:"'
  287. 'CursorDown 1'
  288. 'Put "N:"'
  289. 'CursorDown 1'
  290. 'Put "Pearson r:"'
  291. 'CursorDown 1'
  292. 'Put "r sq.:"'
  293. 'CursorDown 1'
  294. 'Put "Std. Error of Est.:"'
  295. 'CursorDown 1'
  296. 'Put "Intercept (a):"'
  297. 'CursorDown 1'
  298. 'Put "Slope (b):"'
  299. 'CursorDown 2'
  300. 'Put "T test:"'
  301. 'CursorDown 1'
  302. 'Put "Std. Error of Reg. Coef.:"'
  303. 'CursorDown 1'
  304. 'Put "t:"'
  305. 'CursorDown 1'
  306. 'Put "d.f.:"'
  307. 'CursorDown 1'
  308. 'Put "P(T<=t) one-tail:"'
  309. 'CursorDown 1'
  310. 'Put "T-Critical (95%):"'
  311. 'CursorDown 1'
  312. 'Put "T-Critical (99%):"'
  313. 'CursorDown 1'
  314. 'Put "P(T<=t) two-tail:"'
  315. 'CursorDown 1'
  316. 'Put "T-Critical (95%):"'
  317. 'CursorDown 1'
  318. 'Put "T-Critical (99%):"'
  319. 'SelectCell' out_cell
  320. 'CursorDown 1'
  321. 'Column 3'
  322. 'CursorDown 1'
  323. 'Put' NRows
  324. 'CursorDown 1'
  325. 'Put' format(pr,,4)
  326. 'CursorDown 1'
  327. 'Put' format(rs,,4)
  328. 'CursorDown 1'
  329. 'Put' format(sumdifpa,,4)
  330. 'CursorDown 1'
  331. 'Put' format(a,,4)
  332. 'CursorDown 1'
  333. 'Put' format(b,,4)
  334. 'CursorDown 3'
  335. 'Put' format(stdev,,4)
  336. 'CursorDown 1'
  337. 'Put' format(t,,4)
  338. 'CursorDown 1'
  339. 'Put' df
  340. 'CursorDown 1'
  341. 'Put' format(P,,6)
  342. 'CursorDown 1'
  343. 'Put' format(TCRIT1,,4)
  344. 'CursorDown 1'
  345. 'Put' format(TCRIT2,,4)
  346. 'CursorDown 1'
  347. 'Put' format(PT,,6)
  348. 'CursorDown 1'
  349. 'Put' format(TCRIT3,,4)
  350. 'CursorDown 1'
  351. 'Put' format(TCRIT4,,4)
  352. 'CursorDown 2'
  353. 'Column -1'
  354. 'Put "ANOVA:"'
  355. 'CursorDown 1'
  356. 'GetCursorPos'
  357. an_cell=result
  358. 'Put "Source of"'
  359. 'CursorDown 1'
  360. 'Put "Variation"'
  361. 'CursorDown 1'
  362. 'Put "Regression:"'
  363. 'CursorDown 1'
  364. 'Put "Residual:"'
  365. 'CursorDown 1'
  366. 'Put "Total:"'
  367. 'CursorDown 2'
  368. 'Put "P(F Sample <=f) One-tail:"'
  369. 'CursorDown 1'
  370. 'Put "F Critical (95%):"'
  371. 'CursorDown 1'
  372. 'Put "F Critical (99%):"'
  373. 'SelectCell' an_cell
  374. 'Column 1'
  375. 'Put "Sum of"'
  376. 'CursorDown 1'
  377. 'Put "Squares"'
  378. 'CursorDown 1'
  379. 'Put' format(SSReg,,4)
  380. 'CursorDown 1'
  381. 'Put' format(SSRes,,4)
  382. 'CursorDown 1'
  383. 'Put' format(SSTot,,4)
  384. 'CursorDown 2'
  385. 'Put' format(PF,,6)
  386. 'CursorDown 1'
  387. 'Put' format(FCRIT1,,4)
  388. 'CursorDown 1'
  389. 'Put' format(FCRIT2,,4)
  390. 'SelectCell' an_cell
  391. 'Column 2'
  392. 'Put "Mean"'
  393. 'CursorDown 1'
  394. 'Put "Squares"'
  395. 'CursorDown 1'
  396. 'Put' format(MeanSqReg,,4)
  397. 'CursorDown 1'
  398. 'Put' format(MeanSqRes,,4)
  399. 'SelectCell' an_cell
  400. 'Column 3'
  401. 'Put "Degrees of"'
  402. 'CursorDown 1'
  403. 'Put "Freedom"'
  404. 'CursorDown 1'
  405. 'Put' 1
  406. 'CursorDown 1'
  407. 'Put' df
  408. 'CursorDown 1'
  409. 'Put' df+1
  410. 'SelectCell' an_cell
  411. 'Column 4'
  412. 'Put "F Ratio"'
  413. 'CursorDown 2'
  414. 'Put' format(FR,,4)
  415. 'Refresh 1'
  416. 'Refresh 2'
  417. /*'Message' "Finished"*/
  418. /*indicate the main script is finished*/
  419. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  420. result=writeln(6Info, DisplayMsg)
  421. if result~=0 then do
  422.     /*Wait 3 seconds*/
  423.     CALL DELAY(150)
  424.     /* close window*/
  425.     result=close(6Info)
  426. end
  427. 'DEFPUBSCREEN("Workbench")'
  428. exit
  429.  
  430. /* Procedures */
  431.  
  432. cellrow: procedure
  433. do
  434.     parse arg cell
  435.     do charpos=2 to length(cell)
  436.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  437.     end
  438.     return 0
  439. end
  440. Return
  441.  
  442. cellcol: procedure
  443. do
  444.     parse arg cell
  445.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  446.     cell=upper(cell)
  447.     len=length(cell)
  448.     val=0
  449. do charpos=1 to len
  450.     if datatype(substr(cell,charpos,1),n) then
  451.     do cell=reverse(substr(cell,1,charpos-1))
  452.     do x=1 to length(cell)
  453.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  454.     end
  455.     return val
  456.     end
  457.     end
  458.     return 0
  459. end
  460. Return
  461.  
  462. syntax:
  463.      if arg(1)='FAIL' then do
  464.         'Message "Library is unavailable."'
  465.         'DEFPUBSCREEN("Workbench")'
  466.         exit
  467.         end
  468.     'DEFPUBSCREEN("Workbench")'
  469.     exit
  470.  
  471. Format:  procedure
  472.  
  473.      arg number, before, after
  474.      CallLine = SIGL
  475.      if ~datatype(CallLine, 'N') then CallLine = '??'
  476.  
  477.      /* Make sure we have a number as first (required) argument    */
  478.      if ~datatype(number, 'N') then do
  479.         if number = '' then
  480.            fc = 17     /* Wrong number of arguments           */
  481.         else
  482.            fc = 47     /* Arithmetic conversion error             */
  483.         signal FormatSyntaxError
  484.      end
  485.      num = number + 0
  486.      if before = '' & after = '' then
  487.         return num
  488.      else do
  489.         parse var num integer '.' fraction
  490.         if before = '' then before = length(integer)
  491.         if after = '' then after = length(fraction)
  492.         if ~datatype(before, N) | ~datatype(after, N) then
  493.            do fc = 18
  494.            signal FormatSyntaxError
  495.        end
  496.         if before < length(integer) then do
  497.            fc = 18
  498.            signal FormatSyntaxError
  499.         end
  500.         if after ~= length(fraction) then do
  501.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  502.         if integer<1&integer>-1 then integer=integer
  503.            else integer = integer + (fraction % 1)
  504.            fraction = substr(fraction, 3)
  505.         end
  506.         if fraction >= 0 then
  507.            return right(integer, before)'.'fraction
  508.         else
  509.            return right(integer, before)
  510.      end
  511.  
  512.  FormatSyntaxError:
  513.         if show('F', STDERR) then
  514.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  515.         else
  516.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  517.         'Message' mess
  518.         parse source Func .
  519.         if Func = 'FUNCTION' then do
  520.         'DEFPUBSCREEN("Workbench")'
  521.            exit "Err"
  522.         end
  523.         else do
  524.         'DEFPUBSCREEN("Workbench")'
  525.            exit 10
  526.         end
  527.  
  528. PROBT: PROCEDURE
  529.  
  530.     PARSE ARG TA,K1
  531.     A=.36338023
  532.     W=ATAN(TA/SQRT(K1))
  533.     S=SIN(W)
  534.     C=COS(W)
  535.     L=K1-2*INT(K1/2)
  536.     IF L=0 THEN DO
  537.         T1=S
  538.         IF K1=2 THEN DO
  539.         PO=.5*(1+T1)
  540.         RETURN PO
  541.         END
  542.         T2=S
  543.         J1=-1
  544.         J2=0
  545.         K2=(K1-2)/2
  546.     END
  547.     ELSE DO
  548.         T1=W
  549.         IF K1=1 THEN DO
  550.             T1=T1*(1-A*L)
  551.             PO=.5*(1+T1)
  552.             RETURN PO
  553.         END
  554.         T2=S*C
  555.         T1=T1+T2
  556.         IF K1=3 THEN DO
  557.             T1=T1*(1-A*L)
  558.             PO=.5*(1+T1)
  559.             RETURN PO
  560.         END
  561.         J1=0
  562.         J2=1
  563.         K2=(K1-3)/2
  564.     END
  565.     DO I=1 TO K2
  566.         J1=J1+2
  567.         J2=J2+2
  568.         T2=T2*C*C*J1/J2
  569.         T1=T1+T2
  570.     END
  571.     T1=T1*(1-A*L)
  572.     PO=.5*(1+T1)
  573. RETURN PO
  574.  
  575. TCRITICAL: PROCEDURE
  576.  
  577.     PARSE ARG PO,K1
  578.     AO=.0001
  579.     E=.005
  580.     E2=E+E
  581.     A=2*PO-1
  582.     IF K1=1 THEN DO
  583.         TO=TAN(1.57079633*A)
  584.         RETURN TO
  585.     END
  586.     IF K1=2 THEN DO
  587.         SN=SIGN(2/(1-A*A))
  588.         IF SN=-1 THEN TO=-1*(A*SQRT(ABS(2/(1-A*A))))
  589.         ELSE TO=A*SQRT(2/(1-A*A))
  590.         RETURN TO
  591.     END
  592.     P1=PO
  593.     Z=NORMALPP(PO)
  594.     W=Z*(1+2/(1+8*K1))
  595.     T3=K1*(EXP(W*W/K1)-1)
  596.     SELECT
  597.         WHEN Z<0 THEN TT=-1
  598.         WHEN Z=0 THEN TT=0
  599.     OTHERWISE TT=1
  600.     END
  601.     SN=SIGN(T3)
  602.     IF SN=-1 THEN T3=-1*(TT*SQRT(ABS(T3)))
  603.     ELSE T3=TT*SQRT(T3)
  604.     /* LABELA:
  605.     TO=T3
  606.     PO=PROBT(TO,K1)
  607.     F1=PO-P1
  608.     TO=T3+E
  609.     PO=PROBT(TO,K1)
  610.     F2=PO
  611.     TO=T3-E
  612.     PO=PROBT(TO,K1)
  613.     F2=F2-PO
  614.     F2=F2/E2
  615.     T4=T3-F1/F2
  616.     IF ABS(T4-T3)>AO THEN DO
  617.         T3=T4
  618.         SIGNAL 'LABELA'
  619.     END
  620.     */
  621.     T4=0
  622.     DO FOREVER
  623.         TO=T3
  624.         PO=PROBT(TO,K1)
  625.         F1=PO-P1
  626.         TO=T3+E
  627.         PO=PROBT(TO,K1)
  628.         F2=PO
  629.         TO=T3-E
  630.         PO=PROBT(TO,K1)
  631.         F2=F2-PO
  632.         F2=F2/E2
  633.         T4=T3-F1/F2
  634.         IF ABS(T4-T3)>AO THEN T3=T4
  635.         ELSE LEAVE
  636.     END
  637.     TO=T4
  638.     
  639. RETURN TO
  640.  
  641. LOGGAMMA: PROCEDURE
  642.  
  643.     ARG XA
  644.     C1=76.18009173
  645.     C2=-86.50532033
  646.     C3=24.01409822
  647.     C4=-1.231739516
  648.     C5=.001208580
  649.     C6=-.000005364
  650.     C7=2.506628275
  651.     X1=XA-1
  652.     W=X1+5.5
  653.     W=(X1+.5)*LN(W)-W
  654.     S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
  655.     S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
  656.     L=W+LN(C7*S)
  657. RETURN L
  658.  
  659. NORMALPP: PROCEDURE
  660.  
  661.     ARG P0
  662.     A1=2.515517
  663.     A2=.802853
  664.     A3=.010328
  665.     A4=1.432788
  666.     A5=.189269
  667.     A6=.001308
  668.     Q=.5-ABS(P0-.5)
  669.     SN=SIGN(-2*LN(Q))
  670.     IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
  671.     ELSE W=SQRT(-2*LN(Q))
  672.     W1=A1+W*(A2+A3*W)
  673.     W2=1+W*(A4+W*(A5+A6*W))
  674.     ZZ=W-W1/W2
  675.     SELECT
  676.         WHEN (P0-.5)<0 THEN TT=-1
  677.         WHEN (P0-.5)=0 THEN TT=0
  678.         otherwise TT=1
  679.     END
  680.     ZZ=ZZ*TT
  681. RETURN ZZ
  682.  
  683. FCRITICAL: PROCEDURE EXPOSE AC EC EC2
  684.  
  685.     PARSE ARG PO,K1,K2
  686.     /* fIND NORMAL Z CORRESPONDING TO P */
  687.     P1=PO
  688.     Z=NORMALPP(PO)
  689.     IF K2<3 THEN DO
  690.         W=K2**.75
  691.         W=Z/W
  692.         Z=Z*(1.1581-W*(.2296+W*(.0042+.0027*W)))
  693.     END
  694.     /* FIND INITIAL APPROX. TO F */
  695.     A1=2/(9*K1)
  696.     A2=2/(9*K2)
  697.     W=Z*Z
  698.     W1=1+A2*(A2-W-2)
  699.     W2=A1+A2-A1*A2-1
  700.     W3=1+A1*(A1-W-2)
  701.     SN=SIGN(W2*W2-W1*W3)
  702.     IF SN=-1 THEN W3=-1*SQRT(ABS(W2*W2-W1*W3))
  703.     ELSE W3=SQRT(W2*W2-W1*W3)
  704.     FO=(W3-W2)/W1
  705.     FO=FO**3
  706.     /* MODIFIED NEWTON ITERATION TO IMPROVE VALUE OF F */
  707.     DO FOREVER
  708.         FCRIT=FO
  709.         PO=PROBABILITY(FCRIT,K1,K2)
  710.         F1=PO-P1
  711.         FCRIT=FO+EC
  712.         PO=PROBABILITY(FCRIT,K1,K2)
  713.         F2=PO
  714.         FCRIT=FO-EC
  715.         PO=PROBABILITY(FCRIT,K1,K2)
  716.         F2=F2-PO
  717.         F2=F2/EC2
  718.         F1=FO-F1/F2
  719.         IF ABS(F1-FO)>AC THEN FO=F1
  720.         ELSE LEAVE
  721.     END
  722.     FCRIT=F1
  723. RETURN FCRIT
  724.  
  725.  
  726.  
  727. PROBABILITY: PROCEDURE EXPOSE AC EC EC2
  728.  
  729.     PARSE ARG F,K1,K2
  730.     IF K1=1 THEN AC=.00001
  731.     H2=K2/2
  732.     H1=K1/2
  733.     H3=H1+H2
  734.     L1=0
  735.     XX=H1
  736.     LG=LOGGAMMA(XX)
  737.     L1=L1+LG
  738.     XX=H2
  739.     LG=LOGGAMMA(XX)
  740.     L1=L1+LG
  741.     XX=H3
  742.     LG=LOGGAMMA(XX)
  743.     L1=L1-LG
  744.     W1=K2/(K2+K1*F)
  745.     XX=0
  746.     IF H2<(H3*W1) THEN DO
  747.         W2=W1
  748.         W1=1-W1
  749.         W3=H2
  750.         H2=H1
  751.         H1=W3
  752.     END
  753.     ELSE DO
  754.         W2=1-W1
  755.         XX=1
  756.     END
  757.     T1=1
  758.     W3=1
  759.     P=1
  760.     I=INT(H1+W2*H3)
  761.     W4=W1/W2
  762.     /*LABELB:*/
  763.     T2=H1-W3
  764.     IF I=0 THEN W4=W1
  765.     /*LABELC:*/
  766.     DO FOREVER
  767.         T1=T1*T2*W4/(H2+W3)
  768.         P=P+T1
  769.         T2=ABS(T1)
  770.         IF (T2<=AC) & (T2<=AC*P) THEN LEAVE /*SIGNAL 'LABELD'*/
  771.         W3=W3+1
  772.         I=I-1
  773.         IF I>=0 THEN DO /*SIGNAL 'LABELB'*/
  774.             T2=H1-W3
  775.             IF I=0 THEN W4=W1
  776.             ITERATE
  777.         END
  778.         T2=H3
  779.         H3=H3+1
  780.     /*SIGNAL 'LABELC'*/
  781.     END
  782.     /*LABELD:*/
  783.     W1=EXP(H2*LN(W1)+(H1-1)*LN(W2)-L1)
  784.     P=P*W1/H2
  785.     IF XX=0 THEN PO=P
  786.     ELSE PO=1-P
  787.     AC=.001
  788. RETURN PO
  789.  
  790.  
  791.